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Summary 



O 
(N 

Efforts to catalogue the structure of metabolic networks have generated highly detailed, genome-scale 
atlases of biochemical reactions in the cell. Unfortunately, these atlases fall short of capturing the 

.^ kinetic details of metabolic reactions, instead offering only topological information from which to make 

predictions. As a result, studies frequently consider the extent to which the topological structure of a 
metabolic network determines its dynamic behavior, irrespective of kinetic details. Here, we study a 
class of metabolic networks known as non-autocatalytic metabolic cycles, and analytically prove an open 
conjecture regarding the stability of their steady-states. Importantly, our results are invariant to the 

j /^ choice of kinetic parameters, rate laws, equilibrium fluxes, and metabolite concentrations. Unexpectedly, 

our proof exposes an elementary but apparently open problem of locating the roots of a sum of two 
polynomials S = P + Q, when the roots of the summand polynomials P and Q are known. We derive 
two new results named the Stubborn Roots Theorems, which provide sufficient conditions under which 
the roots of S remain qualitatively identical to the roots of P. Our work illustrates how complementary 
feedback, from classical fields like dynamical systems to biology and vice versa, can expose fundamental 
and potentially overlooked questions. 
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"^" Introduction 

O 

Networks of enzyme-catalyzed metabolic reactions are fundamental to the proliferation of life, using the 
energy extracted from environmental nutrients to drive the assembly of organic macromolecules and en- 
able the successful reproduction of the cell. The large-scale architecture of these networks is rich in 
structure: they are broadly organized into overlapping pathways (e.g. catabolic glycolysis and anabolic 
glucoeneogenesis) , and exhibit power-law-like degree distributions with highly connected cofactors (e.g. 



ATP/ADP and NADH/NAD) linking many otherwise-distant metabolites 26 . Perhaps most impor- 
tantly, these networks are capable of robust operation in spite of heterogeneity in the abundances of 
crucial enzymes and substrates [29| . 

To what extent are the robust features of metabolic networks determined by the underlying topolog- 
ical structure of the network itself? This question lies at the center of many studies precisely because 
contemporary metabolic models are largely limited to structural information. Using genomic data, it is 
now possible to reconstruct genome-scale models of metabolism, which predict the presence of absence 
or enzymes (and by virtue, metabolic reactions) in an organism. However, the kinetic details of these 
reactions (such as the rate laws they obey, as well as rate constants like Km or V max ) are largely unknown 
due to the difficulty of measuring them in a high-throughput manner in vivo. As a result, a number of 
generic results linking the qualitative dynamics of a chemical reaction network (such as its potential for 



multist ability or sustained oscillations) to its structural organization 11 12 have begun to populate the 



literature, pointing to a fundamental connection between structure and dynamics. These results make 
remarkably minimal and physically reasonable assumptions on the generic form of kinetic rate laws, that 
render their conclusions largely independent of the choice of parameters. Perhaps the most well-known 



result from studies of this sort is the Deficiency Zero Theorem 11 from chemical reaction network theory 
(CRNT), which gives sufficient conditions for unique equilibria and asymptotic stability of a large class 
of reaction networks. More recent work extending CRNT, such as that of Shinar and Feinberg [29], has 
identified structural properties endowing chemical reaction networks with absolute concentration robust- 
ness (ACR) (i.e. the steady-state concentration of a molecular species is identical in any steady-state the 
dynamical system admits). Importantly, other existing methodologies, including Chemical Organization 
Theory f8ll2l] and the theory of monotone systems [3] , take similar "topological" approaches to under- 
standing how dynamics may be inferred from, and in fact directly influenced by, the structure of reaction 
networks themselves. 



In prior work, we studied a family of non-autocatalytic metabolic cycles 28 , and considered the role 
that their cyclic topology might play in determining their steady-state properties. Our decision to study 
a cycle, as opposed to any other structure, was motivated by the prominent role of cycles (such as the 
TCA and Calvin cycles) in present-day metabolic networks [5] . The approach we took relied on a method 
known as structural kinetic modeling (SKM), which applied a change of variables to the dynamical system 
corresponding to the metabolic cycle. This change of variables enabled us to study the dynamics of the 
metabolic cycle from a structural point of view, with only mild assumptions on the form of the kinetics 
themselves. 



The main outcome of our work in 28 was limited numerical evidence that any steady-state of the 
cycle must be stable to small perturbations, irrespective of equilibrium metabolite concentrations, flux 
magnitudes, or the choice of kinetic parameters. However, we were unable to offer a rigorous analytical 
proof of this claim. In particular, it was unclear whether small regions of parameter space harboring 
unstable equilibria might exist. Perhaps more importantly, computational considerations limited our 
numerical investigations to relatively short metabolic cycles (including up to eight metabolites), leaving 
open the possibility that instability appeared as cycles grew longer. As a result, we left the question of 
stability as an unproven conjecture (herein referred to as the cycle stability conjecture). This conjecture 
is the object of study in the first part of this work. 

The difficulty with proving the cycle stability conjecture reduced to locating the roots of a high- 
order polynomial in the complex plane. Although a number of classical results from control theory are 
commonly applied to problems like this, they were rendered largely unusable for the polynomial in [28] . 
For example, the Routh-Hurwitz (RH) criterion, perhaps the best-known technique for constraining the 
locations of a polynomial's roots, requires the precise calculation of coefficients of the polynomial under 



study. Unfortunately, the calculation of these coefficients for the polynomial in 28 became analyti- 
cally intractable as the number of reactions in the cycle (which we assumed to be arbitrarily large) 
grew. Furthermore, because these coefficients were themselves functions of SKM variables, and were 
only constrained to lie in complicated intervals on the real axis, the problem of proving stability became 
substantially more difficult. Well-known methods, such as Kharitonov's theorems [7J, exist to study how 
uncertainty in the coefficients of a polynomial impacts its roots. In addition, the Method of Resultants, 



used by Gross and Feudel in the context of studying generalized models 19 , can be used to identify 



when pairs of imaginary roots cross the imaginary axis (suggesting the onset of instability and sustained 



oscillations). However, the complexity and scale of the polynomial in 28 again rendered the application 
of both Kharitonov's theorems and the Method of Resultants infeasible. 

Our difficulty with bringing classical tools to bear on the cycle stability conjecture motivated us 
to revisit the problem from a completely new perspective. In this work, we resolve the cycle stability 
conjecture by reformulating it as a question of locating the roots of a sum of two polynomials. By doing 
so, we are able to apply a classical technique (Rouche's Theorem [6]) from complex analysis to resolve the 
conjecture and prove the stability of non-autocatalytic metabolic cycles. Quite unexpectedly, our proof 



leads us to a substantially more general question: how do the roots of a sum of polynomials S = P + Q 
depend on the roots of P and Q themselves? Using a method identical to the one used to prove the cycle 
stability conjecture, we prove two new theorems, which we call the Stubborn Roots Theorems, which 
give sufficient conditions for when the roots of S are qualitatively identical to the roots of P. To our 
knowledge, there are few generic results which provide information regarding the locations of the roots 
of a sum of two polynomials. Given the fundamental importance of locating the roots of a polynomial 
in the study of dynamical systems (and in applications of dynamical systems to fields such as in systems 
biology), we feel that the Stubborn Roots Theorems may find use in other contexts. 

Results 

Stability of Metabolic Cycles 

We begin by presenting a model of the dynamics of a simple non-autocatalytic metabolic cycle, and then 
proceed to using SKM to study the stability of its equilibria. First, we describe the generic structure of the 



metabolic cycle under study, which is identical to the one studied in 28 . The cycle contains n metabolites 
(Mi . . . M n ) and two cofactors (0\ and 2 ) which provide the energetic force to thermodynamically drive 
the metabolic reactions, and can be illustrated by 

— > Mi 
Mi + Oi — > M 2 + 2 

Mi — > M i+ i, i = 2...n- 1 

M„ — ► Mi 
M n — ► 

o 2 ^^ 1 (1) 

In this cycle, each metabolite Mi is converted to metabolite Mi+i for i = 1 . . . n — 1. A constant flux 
of Mi enters the system. A proportion of the last metabolite M„ is converted back to Mi, while the 
remainder leaves the system. The high-energy cofactor 0\ is converted to its low-energy cofactor partner 
2 in the reaction catalyzing the conversion of Mi to M 2 . In a separate reaction, energy is input into 
the system to drive the reformation of the higher energy molecule 0\ . 

At steady-state, the magnitude of the flux through each reaction in the cycle can be calculated by 
enforcing mass balances on each metabolite. To do so, we assume that a constant flux of generic magnitude 
av,0 < a < 1 of metabolite Mi flows into the network. A proportion (1 — a)v of the flux entering M„ 
is channeled back towards M 1; while the remaining flux av exits the system. All other reactions carry 
a steady-state flux of v. It is easily verified that this flux vector is in the nullspace of the stoichiometric 
matrix S (see Appendix). We assume that the kinetics of each reaction are monotonic, that is, that 
an increase in the concentration of substrate for any reaction will consequently increase the rate of the 
reaction. This assumption is quite generic, and is amenable with many well-known biochemical reaction 
mechanisms, including the law of mass-action as well as Michaelis-Menten and Hill kinetics. 

To prove the stability of a steady-state of ([!]), we must prove that the Jacobian of (Tij), evaluated at an 
arbitrary steady-state, always has eigenvalues with negative real part. As shown in |28| (and re-derived 
in the Appendix, Equations (25|27)), the Jacobian J n for the metabolic cycle of size n illustrated above 



can be calculated using SKM to be 



■J n 
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7 n+2 — "n+3- 



(2) 



Crucially, the assumption of monotonic kinetics constrains all the elasticities, fi'i, in (|2|) to be greater 
than zero (see Appendix). In 28 , we provided evidence J n could not have eigenvalues with positive 
real part. Below, we proceed to analytically prove this conjecture. To simplify some calculations, we 
elect to work with the negative counterpart of the Jacobian, J~ = — J n , and prove that J~ cannot have 
eigenvalues with negative real part. This is equivalent to proving that J n cannot have eigenvalues with 
positive real part. 

First, we calculate the characteristic polynomial of J~ , which we call Xn(A), explicitly (calculations 
shown in Appendix): 



P n = ((0i- X)(9 n +2 + o n+3 - A) - 9 1 6 n+2 ) (0„ + e n+1 - A) J] ($i - A) 

i=2...n-l 

Qi=-(e n+3 -\) n °i 

i—l...n 
Xn(\)=Pn+Ql 



(3) 

(4) 

(5) 



Thus, Xn is the sum of two polynomials, an n + 1 th order polynomial P n and a first-order polynomial 
Qi (note that the subscript n denotes the size of the cycle, not the degree of the polynomial). Next, we 
prove three lemmas on the relative location of the roots of P n and Qi, showing that they are strongly 
constrained. Later on, these constraints will be crucial to proving that Xn can only have roots with 
positive real part. 

Lemma 1. All of the roots of P n are positive and real. 

Proof. By inspection, at least n— 1 of P n 's roots, contained in the product term of (pj), must be positive 
and real. For the remaining two roots, we must study the quadratic polynomial 



(0i -A) (ft 



n+2 



%i+3 — A) — 



'lCn+2- 



(0) 



First, we prove that the roots of ([6]) must be real. Calculating the discriminant A of this quadratic 
polynomial, we find 



A — (01 + 9 n+ 2 + 0n+3) ~ 40i0 n+ 3 



J U ' 0„+3 + 20ift„ + 2 + 20i0 n+3 + 20 n+2 0n+3 — 40i0 n +3 
•i- 1 •> 02 , 2 



1 T w n+2 T u n+3 

I — 28 1 8 n+ 3 + n+3 + 8 n+2 + 20 n+2 n +3 + 20i0„ +2 
2 , n2 



— (01 — dn+3) + S n +2 + 20 n +20n+3 + 20i0„ + 2 

>0. 
Since A > 0, ffify cannot have imaginary roots and all of the roots of P n are purely real. 



Next, we show that the pair of roots of ffib must be positive. If we expand (|6]), we find 

a 2 - A(0i + e n+2 + e n+3 ) + Mn+3 - o. (7) 

The product of the two roots of Q are 9i9 n +3 > 0, and the sum of the roots is 6\ + 9 n+ 2 + d n +3 > 0. 
Therefore, both of the roots of ffh must both be positive. 

□ 

We next prove a related lemma regarding the location of the root of Q±. 

Lemma 2. The root r q of Qi must be larger than at least one root of P n , and smaller than another root 

0fP n . 

Proof. Consider the quadratic factor of P n , {9\ — \)(0 n +2 + #n+3 — A) — 6i& n +2- By Lemma 1, this 
quadratic polynomial has distinct (since A > 0) real, positive roots. Let pi,p 2 denote these roots, 
ordered by magnitude so that p\ < p 2 - Since the leading term of the quadratic is positive, the roots 
of the polynomial divide the real line into 3 regions: {A < p{\ and {A > p 2 } where the polynomial is 
greater than or equal to zero, and {p\ < A < p 2 } where the polynomial is strictly negative. By inspection, 
f q = n +3- Directly evaluating the value of the quadratic polynomial at A = r q — #„+3, we find: 

®n+3 ~ 6n+s{9l + Qn+2 + ^n+3j + #l#n+3 = _ &n+3&n+2 < 0. (8) 

Because the value of the quadratic factor is negative at A = r q , r q must lie between the roots of (m). 

a 

Finally, we prove a lemma regarding the magnitude of P n and Q\ at the origin. 

Lemma 3. |P„(0)| > |Qi(0)| 

Proof. Explicitly calculating P n (0), we find |P n (0)| = \(8 n + 9 n +x)6 n +s IXi=i n-i^'l- This is always 
greater than |Qi(0)| = |^n^n+3lli=i n-i ^1- 

n 

Now, using Lemmas 1-3, we proceed to prove that the roots of x-n must lie in the positive real half of 
the complex plane. To do so, we will make use of a well-known theorem from complex analysis known as 
Rouche's Theorem. 

Theorem 1 (Symmetric Rouche's Theorem). Two holomorphic functions f and g have the same number 
of roots within a region bounded by some continuous closed contour C (on which neither f nor g have 
any poles or zeros) if the strict inequality 

\f(z)-g(z)\<\g(z)\ 
holds on C fm. 

In essence, Rouche's theorem offers a way to determine the number of roots of a difference of two 
functions lying inside a closed contour in the complex plane. We will use / = S = P n + Q\ and g = P n , 
By taking contours bounding larger and larger regions of the left half-plane (those complex numbers with 
negative real-part), we will show that \Qi\ < \P n \ on the contour, and thus prove that S has no roots 
inside this contour, i.e. no roots in the left half-plane. 

We let our contour C^ consist of two parts (see Figure 1): 



• the portion of the circle {|A| = R} centered at the origin in the negative real half of the complex 
plane 

• the portion of the imaginary axis connecting the two points of intersection of the above circle with 
the imaginary axis 

IP I 
First, we will prove that along an arc of sufficiently large radius, j^pj > 1. Since n > 0, given an 

arc of sufficiently large radius R, \Qi\ < \P n \ simply because P n is of higher order than Q 1 (i.e. the 

highest-order term in P n is A" +1 which dominates the highest order term in Q 1; A 1 , for very large A). 

ip i 
Next, we will prove that along the upper half of the imaginary axis, jtjH > 1. To do so, let us consider 



IQil 
aking 
of P n as r, for i = 1, . . . , n + 1 and the root of Qi as i" q ,*we have 



i p i 
the behavior of 1^4 by substituting A = iy, y > into (|5|) and taking the modulus. Denoting the roots 



\ p n(iy)\ = \(n - iy){r 2 -iy)... (r n +i - iy)\ 

= ^{ri + V 2 ){rl + y 2 )...{rl +1+ if) (9) 

\Qi(iy)\ = \c-iby\ 

= ^c 2 + b 2 y 2 = b^r 2 q +y 2 , (10) 

where c = |Qi(0)| and r q = f . We must show that the following condition holds for all y > 0: 

IP I J{rl+y 2 ){rl+y 2 )---{rl +l +y 2 ) 

x{y) = \^\ = * r == > 1. (ii) 

\Qi\ bjrl + y 2 



Note that at y = 0, we know (11 1 is satisfied because |P„(0)| > |Qi(0)|. If we can show that x(y) is a 
strictly increasing function, then we know that pTj ) will be satisfied for all y. To do so, let us 
x 2 (y). Note that x(y) > 0, x 2 (0) > 1, and if j-ar> for all y, then x 2 (y) > 1 for all y. This i 



dy 

imply that x(y) > 1 for all y > 0. We have 



work with 
would then 



2/ x _ \P^_ _ (r 2 i+y 2 )(r 2 2 +y 2 )...(r 2 n+1 +y 2 ) 

Taking a derivative of x 2 with respect to y and using the identity ^(/l/g) = ~ d ^f2 + -^ 2 -fi, where 
A = \Pn\ 2 ih = iqti (the product rule applied to \P n \ 2 and l/|<5i | 2 ), we find 

u 1 y \i=l...n+l \jjLi ) 1 y fe=l...n+l 

Without loss of generality, suppose r\ is the smallest root of P n . From the first term on the right 



hand side of (13), select the term corresponding to i = 1. Recall that, by Lemma 2, r\ must be smaller 
than r q . Isolating just this term and the negative term in the parentheses on the right-hand-side of (13) 
and summing, we find 



r% + y 2 



1/2 ] n ^ + y 2 - ( 14 ) 



1 y 7 i=2...n+l 



Since r q > ri, (14 1 is positive. There are no more negative terms in (13), proving that x (y) is 
strictly increasing. This proves that \P n \ > \Qi\ on the positive imaginary axis. Furthermore, since real 
polynomials are symmetric across the real axis, an identical argument shows that \P n \ > \Qi\ on the 
negative imaginary axis. In particular, setting A = — iy for y > yields the exact same expressions for 
\P n \ and |Q X |. 

We have satisfied all of the assumptions of Rouche's Theorem, and have proven x n (A) contains no 
roots in the left half of the complex plane. Therefore, the nonautocatalytic metabolic cycle always has 
stable equilibria. 

The Stubborn Roots Theorem 

Can we use the methods illustrated in the prior section to locate the roots of the sum of two more 
general polynomials? Our motivation for studying this problem derives from control theory, where it 
is common to ask whether the roots of a polynomial lie in one half of the complex plane Hi. Such 
polynomials frequently correspond to the characteristic equation of the Jacobian matrix of a dynamical 
system. Although we do not provide a generic method for predicting whether a matrix's characteristic 
equation may be written as the sum of two simpler polynomials, the appearance of such structure in our 
studies of a metabolic cycle suggests that related, "well-ordered" systems may exhibit similar properties. 
The main question we ask in this section is under what conditions may the roots of a polynomial P 
be "stubborn" when P is summed with another polynomial Q: in such a case, the roots of the summed 
polynomial S — P + Q remain qualitatively identical to those of P. By qualitatively identical, we mean 
specifically that the number of roots of P in the left (right) half of the complex plane is equal to the 
number of roots of S in the left (right) half of the complex plane. This question follows in the spirit 
of similar work by Anderson 2 . Our primary result is a theorem, which we call the Stubborn Roots 
Theorem, which gives sufficient conditions under which the location of the roots of a sum of polynomials 
S = P + Q remains qualitatively unchanged from P. 

Theorem 2 (Stubborn Real Roots Theorem). Let P n and Q m be polynomials of order n and m, respec- 
tively, andletn>m. Assume that all the roots of P n andQ m are purely real, and that |P„(0)| > |Q m (0)|. 
Denote by Pi and qi the roots of P n and Q m ordered by magnitude, so that \pi\ < \p 2 \ < ... < \p n \ and 
I Si I < \l2\ < • • . < \q m \- If for every j = l...m, \pj\ < \qj\, then the number of roots of S — P n + Q m 
located in the negative (positive) real half of the complex plane is equal to the number of roots of P n in 
the negative (positive) real half of the complex plane. 

The proof of Theorem 2 is provided in the Appendix, and follows precisely the same line of reasoning 
as the proof of the cycle stability conjecture in the previous section. The theorem relies on two critical 
assumptions relating P n and Q m . First, at the origin, |P n (0)| > |Q m (0)|. Second, it must be possible to 
assign to each root of qi of Q m a unique root pi of P n such that qi > p%. 

The power of Theorem 2 is that it enables one to qualitatively locate the roots of a polynomial S 
simply by inspecting the roots of its summands P n and Q m . If the roots of P n and Q m are easily 
calculated (as in the case of the cycle stability conjecture in the prior section), then the roots of S can be 
immediately located without resorting to difficult calculations. In many ways, Theorem 2 is reminiscent 



of the work reported in 14 . There, Fisk describes the behavior of the roots of sums of polynomials which 
"interlace." For two polynomials P and Q to interlace, the roots of P and Q alternate when ordered from 
most negative to most positive, so that p\ < q% < p 2 < qi ■ ■ ■■ Notably, our result here is more general, 
and includes interlacing as a special case. 



Stubborn Complex Roots 

What happens when matters become complex? In this section, we generalize the Stubborn Roots Theorem 
to cases when the roots of P are not necessarily all real. This is often the case in dynamical systems, 
where complex roots indicate oscillatory phenomena such as spiraling or limit cycles [31]. Proceeding 
along the same lines as before, we find that the roots of P remain stubborn to the addition of Q as long 
as they remain predominantly real. That is, if the real component of the complex roots of P is larger 
than their imaginary component, then a more general version of the Stubborn Roots Theorem holds. 

Theorem 3 (Stubborn Complex Roots Theorem). Let P n and Q m be polynomials of order n and m, 
respectively, n > m. Let the m roots of Q m be positive and purely real. Further, let P n have at least m 
real roots, and let the remainder of the roots be either real or complex. Assume that \P n (0)\ > |Q m (0)|. 
Furthermore, assume that for each complex root of P, pk, the magnitude of the real component \Re(pk)\ 
is larger than the magnitude of its imaginary component, \Im(pk)\- Denote by pi and qi the real roots of 
P n and Q m ordered by magnitude, so that \pi\ < \p 2 \ < ... < \p m \ and \Qi\ < 1 92 1 < •■■ < \lm\- If for 
every j = l...m, \pj\ < \qj\, then the number of roots of S = P n + Q m located in the negative (positive) 
real half of the complex plane is equal to the number of roots of P n in the negative (positive) real half of 
the complex plane (Figure 2). 

Proof. We prove the theorem for the case when P n has n — 2 real roots, two complex conjugate roots 
p+ = a + bi and p~ = a — bi, and Q has one positive real root. The result can be straightforwardly 
(via wrenching and tedious algebraic calculations) extended to the generic case in Theorem 3 using an 
identical argument. As before, we apply Rouche's Theorem using a half circle in the negative real half 
of the complex plane using f = P = P n + Qi and g — P n . First, we consider the behavior of the two 
polynomials on the large arc in the negative real half of the complex plane. As before, \P n \ dominates 
|Qi| as the radius of the arc grows larger. 

Turning our attention to the behavior of the polynomials on the positive imaginary axis, we substitute 
z = iy to find 



\Pn\ 2 = U7^(y 2 +pf)((y - b) 2 + a 2 )((y + bf + « 2 ) 
IQil 2 y 2 + ql 



rbrii = lu=i l ' \,\ 2 — ' (15) 



where for simplicity of notation, we have assumed P n has leading coefficient 1. Differentiating and 
simplifying algebraically, we obtain the expression 



+ (2(y-b){(y + bf + a 2 ) + 2(y + b){{y-bf + a 2 )) n ~fl{p 2 + y 2 ), (16) 

where C — ((y — b) 2 + a 2 )((y + b) 2 + a 2 ). Using an argument identical to the one used to prove the 



cycle stability conjecture, it is clear that the top term in ( 16 1 is positive. We are then left to ensure that 



the bottom term is positive. Expanding this term, we find 

n-2 



~ 2 ■ V) 



(2{y-b){(y + b) 2 + o 2 )+2{y + b)({y-b) 2 +a 2 ))\{{p 2 

i=l 
n-1 

= (4y(y 2 + a 2 + b 2 ) + 2b(-4by)) J[ (p 2 + y 2 ) 

i=\ 
n-2 

= (4y(y 2 + a 2 ~b 2 ))Y[(p 2 +y 2 ). (17) 



Thus, (17) is certain to be positive as long as \a\ > \b\. In this case, we can once more apply Rouche's 
theorem to show that S = P n + Q\ has the same number of roots in each half of the complex plane as 

n 
Discussion 

A major challenge in systems biology is efficiently studying biological networks through detailed atlases 
of their topological structures. These atlases, assembled from the cumulative results of many high- 
throughput, large-scale experiments, capture many of the physical links which underlie such networks. 
However, they fail to describe most of the detailed dynamics taking place on the network itself. Here, we 
have studied the stability properties of a generic type of metabolic network, and analytically proven that 
under quite mild assumptions on the reaction kinetics, any steady-state of the network must be stable. 
To prove this, we re-formulated the question of stability as a problem of locating the roots of a sum of two 
polynomials whose roots were easily calculated. This reformulation exposed a more fundamental problem 
of locating the roots of the sum of two polynomials, and we proved two new results (the Stubborn Roots 
Theorems) offering sufficient conditions under which the roots of the sum are not qualitatively different 
from the roots of one of the summands. 

The study of metabolic networks and their stability has played an important role in research into the 
origin of life. Many of the earliest papers studying simple models of primordial metabolic networks fo- 
cused on elucidating their stability properties, hypothesizing that molecular self-organization may arisen 
through self-sustaining metabolic cycles [9]. Because early metabolism almost certainly lacked the com- 
plex regulatory mechanisms and circuits which appear in cells today, these investigations suggested that 
stability of these cycles to fluctuations in environmental conditions was necessary for their survival. Pur- 
suing this line of thought, Piedrahta et al recently proposed a simple chemical reaction network composed 
of interlocking cycles which could establish and maintain a stable steady-state, even in the face of a sud- 



den loss of some constituent metabolite of the cycle itself 25 . In response to such a catastrophe, the 
remaining metabolites re-produced the missing metabolite. This notion of closure, in which all of the 
metabolites necessary for sustaining the system can be produced by the metabolic network itself, also 
plays an important role in Chemical Organization Theory, a method distinct from SKM for studying 
chemical reaction networks based on structure and mentioned earlier |8| |21) . 

Naturally, one may ask: how important is the stability of equilibria to the robust function of present- 
day biological networks? A rich and diverse literature, dating back nearly half a century and still ex- 



panding today, describes the importance of stabilizing structures in ecological networks 20p4 . However 



the extent to which stability plays a role in the fitness of metabolic systems is unclear, and studies inves- 
tigating this question have failed to produce a definitive conclusion. In particular, if stability endowed 
metabolic networks with some evolutionary advantage, one should expect an enrichment for stabilizing 
features and structures in contemporary metabolic networks. While some studies have demonstrated the 



10 



importance of some key stabilizing edges in metabolic networks, such as the allosteric feedback of ATP 



onto phosphofructokinase in glycolysis 17 30 , others have failed to identify an enrichment of stabilizing 
structures in metabolism as a whole [32). In fact, synthetic biologists routinely exploit instability in 



order to generate circuits exhibiting sustained oscillations, both in metabolic 16 and transcriptional 10 
systems. Importantly, we note that the work presented here only studied dynamics near equilibrium 
points, and ignored nonlocal dynamics (such as the appearance of periodic orbits arising from global 
bifurcations). Efforts extending generalized modeling and SKM to understanding nonlocal dynamics, are 



now appearing in the literature 22 



We expect that the results presented here may find useful application in several challenges facing 
contemporary biology. First, SKM and generalized modeling could be used as coarse-grained techniques 
for vetting synthetic circuit designs for their potential to exhibit desirable behaviors (such as robust 
stability). In a prior study [27] , we did precisely this, using generalized modeling to identify which 
topological circuit designs were entirely incapable of oscillations, irrespective of the choice of kinetic 
parameters or rate laws. Second, a great deal of interest now exists in using high-throughput metabolomics 
and fluxomics data to identify the role that small-molecule (e.g. allosteric) regulation of metabolic 
enzymes plays in shaping metabolic dynamics [23] . Given the difficulty in accurately measuring kinetic 
rate constants in vivo, we envision that SKM (and our results here highlighting the inherent stability 
of certain topological motifs) might serve as a useful bridge between detailed mechanistic models and 
experimental data, highlighting those regulatory interactions which are crucial to the robust function of 
the network as a whole. 

Finally, our results (in particular, the Stubborn Roots Theorems) illustrate the potential for biological 
questions to reveal interesting and unsolved problems in other fields. What appeared to us initially as 
a simple problem of locating the roots of the sum of two stable polynomials, quickly blossomed into the 
exploration of widely diverse fields of active research, from control theory to matrix analysis. There is 
now a growing number of examples of similar feedback from biology to other fields, from classic results 
in evolutionary optimization [15 to to the design of novel algorithms [I]. Interestingly, many of these 
cross-fertilizations of ideas have taken place because of an abundance of biological data, and a need for 
analytical tools to understand it. Here, it has been quite the contrary: our study of a topological model of 
a metabolic cycle was motivated by a dearth of data on the kinetics of metabolic reactions. Nevertheless, 
in both cases, the ultimate outcome is deeper understanding, relevant to both biology and the fields from 
which it draws new tools and ideas. 

Appendix 

Structural Kinetic Modeling 

To study the stability of a steady-state of a metabolic network, we employ a technique known as structural 
kinetic modeling (SKM) [ 30] . SKM is a non-dimensionalization procedure which replaces conventional 
kinetic parameters (such as V max and Km) with normalized parameters known as elasticities. In the past, 
SKM (and its generalization, known as Generalized Modeling (GM)) has been paired with complementary 
methods studying other dynamic features of a system, such as the effects of noise [17] . 

As illustrated below, elasticities have several properties which make them powerful tools for studying 
metabolic dynamics. First, in contrast with kinetic parameters (whose values may be uncertain over many 
orders of magnitude), elasticities are constrained to lie in well-defined ranges (for example, between zero 
and one), and sampling elasticities across this range effectively captures all possible values of kinetic 
parameters. Second, the value of an elasticity does not depend on the particular choice of kinetic rate 
law. Instead, an elasticity is simply a normalized measure of the sensitivity of a rate law to infinitesimal 
changes in a metabolite's concentration. 

To study the stability of a steady state, SKM calculates the Jacobian matrix J of the dynamical 
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system corresponding to a metabolic network. If we let C be the m-dimensional vector of metabolite 
concentrations, N be the m x r stoichiometric network, and v be the r-dimensional vector of metabolic 
fluxes, then the dynamics of a metabolic network are governed by the system of differential equations 

^ = Nv(C,k), (18) 

where k is a vector of parameters and v(C,k) indicates that the vector of fluxes is dependent on both 
metabolite concentrations and kinetic parameters. Assuming that a non-zero steady-state C° exists, we 
can make a change of variables and write 

X -^- A-=N- V ^ C °^ »-(x)- ^ (C ' k) (19) 

' ~ C° 3 ~ J C? AM X 7 — V .(C° k)' 

where i = 1 . . . m and j = 1 . . . r. 

Then, we can write the Jacobian as 

J = A 9 ^ = A0. (20) 

The stability of the steady-state C° is then dependent on the eigenvalues of J. If the real component of 
all eigenvalues of J are negative, then the steady-state is stable. Thus, the problem of stability reduces 
to finding the eigenvalues of J. 

The element which encodes the effective kinetic dependence of reaction rates on metabolites of SKM 
is the r x m elasticity matrix 0. The (i, j) element of describes the sensitivity of the normalized rate 
of reaction i to the normalized concentration of metabolite j. This corresponds precisely to the effective 
kinetic order of the i th reaction with respect to the j th substrate: if the rate of reaction is linear with 



the amount of substrate, then 6 = 1, while if it is zeroth order, 9 = 13 . Importantly, we assume that 
all elasticities in the metabolic cycle are greater than zero; that is, that an increase in the substrate of 
any reaction will increase the rate of that reaction. The analytical power of SKM comes precisely from 
the constrained and well-defined ranges of each element of . 

To illustrate the utility of elasticities, we derive below the elasticity of a metabolite involved in a 
Michaelis-Menten reaction. Consider a biochemical reaction governed by the rate law 



K M + S 
Assuming that this reaction is embedded within a reaction network where metabolite S is at equi 



librium concentration So, we can calculate the normalized reaction rate \x by normalizing (21 1 by its 
steady-state reaction rate: 



-v max s 

K M +S 
Km+S 

S(K M + go) 

S (K M + S) 

x(K M + go) 

K M + a;<So 
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where x — S/ So is the normalized concentration of S. Then, the elasticity is 



()fl: 

d.r 



1 



K M 



x=l 



K 



M 



Sn 



(23) 



Notice that since Sq > 0, 6 is constrained to the range (0,1). The outcome of applying SKM to 
an entire metabolic network is a Jacobian, whose elements are formulated in terms of elasticities with 
well-defined ranges. Prior studies have used computational surveys |18[|30| as well as analytical work 28 
to study the role that particular key elasticities play in determining the stability of the network. 

Finally, it may be useful to give a bit more intuition regarding the generality of an elasticity. To do 
so, we consider a reaction, governed by Michaelis-Menten kinetics, which exhibits an elasticity (explicitly 
calculated in Equation p3| above) of 0.5. First, note that this elasticity may correspond to any combi- 
nation of So and Km which satisfy -g — ¥g- =0.5, for example Km = So = 1 or Km = So = 2. Thus, a 
single value for an elasticity in fact corresponds to a large locus of steady-state concentrations So- Now, 
notice that if we consider all elasticities in the range (0, 1), we in fact capture all possible combinations 
of So > and Km > 0! Thus, if we prove a theorem using this elasticity, and this theorem holds for all 
values of the elasticity in the range (0, 1), then it similarly holds for all possible choices of So and Km- 
Obviously, the result also holds for any other kinetic rate laws for which this elasticity is valid. Thus, 
we effectively capture the entire space of possible parameters and steady-state concentrations. This is 
precisely the approach taken in proving the cycle stability conjecture. 

Characteristic Polynomial for the General Metabolic Cycle 

Here, we formulate the structural kinetic model for the metabolic cycle depicted in (IT]) and reproduced 
below. 



— > Mi 
Mi + Ox — > M 2 + 2 

M % —fM l+ll i = 2...n-l 
M„ — > Mj 
M„ — * 

energy 
Ui > U\ 



(24) 



The analysis presented below is identical to that presented in 28 . The stoichiometric matrix S of 
the metabolic cycle is 



S = 



1 -1 



1 







1 -1 ••• 
1 •■• 





0-10 
1 



-10 
1-1-10 
1 
0-1 



(25) 
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where the rows correspond to each metabolite in the system. The generalized forms of A, the normalized 
stoichiometric matrix, and ©, the elasticity matrix, are shown below for the system depicted in (nl). Note 
that we make no assumptions on the steady state concentrations of metabolites, denoted by the vector 
(Mi,M2, . . . ,Mjv, 0\, O-i). Furthermore, the steady-state flux through the cycle is equal to a generic 
magnitude v, v > 0, except for the input and outflow reactions (with flux av) and the reaction from M n 
to Mi (with flux (1 - a)v). 



a i ■ 


— V 
Mi 


•• 








U-o) 





Mi 


Mi 





V 

Mi 


— v 
M 2 




















V 















and 










— V 



Mjv-i 

V 

M N 







— av 
M N 







(1 — Oi)v 

M N 





(26) 







On 

^N+l 






^N+2 







9jv+3Qi 

o 2 



(27) 



With N + 1 metabolites and N + 3 reactions, A is N + 1 x (N + 3) and © is (N + 3) X N + 1. Note that 
the last row (corresponding to cofactor O2) of A is omitted because the cofactors come as a conserved 
pair, and the bottom right element of (corresponding to the dependence of the last reaction on O2) 
is replaced with a negative element in order to a cco unt for this conservation (for more information on 



modeling of conserved moeties, see the SI Text of 

the eigenvalues of the negative of J n , which we call J, 



30 



Then, the Jacobian J n = A0. We elect to study 
". Note that the eigenvalues of J r 7 are precisely 
the negative of the eigenvalues of J n . Thus, proving that all the eigenvalues of J~ have positive real 
part is equivalent to proving that all of the eigenvalues of J n have negative real part. The characteristic 
polynomial of J~ , x„(A) is 



Xn(A) 



M 



- -X 










-81 

M 2 


02 

M 2 


A . 








-6 2 
M 3 




0,1-1 

• AT„_i 










-0«-i 
M„ 


0i 
Ox 











— On 
Ml 









e n +e n+ i 
m„ 



-A 



Mi 

M 2 



2+0 „ 



Ox 



Above, we have made a change of variables so that 6i = 6i,i — 2...71+2, B\ = o&\, and 9 n +3 = — 
Note that this change of variables does not affect our assumption that all 9 are greater than zero. 
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Now, we proceed to explicitly calculate Xn(A) by calculating the determinant. Expanding along the 
n — 1 th column, we have: 



Xn(A) = (- 







Si 

Mi 

-0i 

M 2 


A 

312- 7 A 








Mi 






^77 + 2 

Mi 

— #n + 2 

M 2 


#71-1 


A) 





Oi 


-02 

M 3 





e n -2 

•■• M n _ 2 





A 




M„ A 







#7i + 2 + #77 + 3 

Oi 




Th ~ A 
-Si 
M 2 




e 2 \ 

372-~ A - 






~6„ 

Mi 




071 + 2 

Mi 

— 0ti + 2 

M 2 






On-1 

M n 





-e 2 

M 3 


On -2 \ 

M„_ 2 A 


















Oi 






-e>^-2 

Af n _i 










$n + 2+0?i + 3 

Oi 


-A 





The important thing to note is that now, the n — 1 th column in the first matrix and the n th row in 
the second are zero except for a single entry. We can continue expanding along such columns until we 
arrive at 



X»(A) = (^-A)...(Tr-A) 



Mi 



A 





lx 
Oi 



'71 -1 



M„ M 3 



Mi 



M 2 



Mi 





Mi 

"m; a 



^77 + 2 

Mi 

— #77 + 2 
M 2 

-1 + 2+^77 + 3 



o 



"ti + 2 

"mT 


i + 2+0n + 3 

Oi 



Finally, we have, for n > 3: 



XnW = I (]y--A)( A) 



Oi 



( 



'n+3 

Oi 



a ) n 



A/, 



Mid 



'lCn+2 \ yfn T Cn + 1 



M„ 



A ) II (if-A) (28) 



i=2...n-l 



If we make a final change of variables so that 0, = ^r, j = 1 . . . n, #„+i = ^ 1 , #„+2 = q +2 , # n +3 
— , then we arrively precisely at &5 



Characteristic Polynomial for a Simple Metabolic Cycle 

To give the reader more intuition regarding the mechanics of SKM calculations, we briefly describe 
how SKM may be used to model the dynamics of a simple, nonautocatalytic metabolic cycle with two 
metabolites and one cofactor pair. The reactions of this cycle are 
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0^ A 
A + Ox—tB + 02 

B — ► A 
5^0 

„ energy _ 
C2 ^ L^l 



(29) 



Then, the stoichiometric matrix S, the normalized stoichiometric matrix A, and the elasticity matrix 
©, may be written 

1-11 
1-1-10 
0-10 1 
1 0-1 



(30) 



A = v 
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-1 

A" 




A" 













1 
B° 
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— a 

B<> 







B° 







-1 
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1 
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9i 







^4 






= 






o 2 

0:; 










' 














9 5 . 





(31) 



(32) 



where v is an arbitrary unit of flux, a £ (0, 1), and ^4°,i3 , and 0° are the steady-state concentrations 
of A, B, and 0\, respectively. 

Note that A has precisely the same structure as S, but with one fewer row (due to the mass conser- 
vation associated with the cofactor pair, discussed in the prior section). Finally, the Jacobian for this 
system may be straightforwardly written 



A© = v 



-01 

A" 


(l-a)0 2 

A° 




-0. 

A" 


0i 

B" 


-(l-a)« 2 - 
B° 


n:0 :i 


04 

BO 


0i 







-04 

Of 



(33) 
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Proof of the Real Stubborn Roots Theorem 

Consider the sum of two polynomials 

S = P n + Q m (34) 

P n = (Pi - A)(P2 - A) . . . ( Pn - A), n > 2 (35) 

Qm = (gi - A)(g 2 - A)...(g m - A),m < n. (36) 

Above, all roots pi and g, are real and positive. Further, we assume that |P„(0)| > |Qi(0)|. We claim 
that P cannot have any roots lying in the left half of the complex plane. To prove this, we make use of 
the symmetric form of Rouche's Theorem (Theorem 1 in the text) We will use / = P = P n + Q m and 
g = P„. Our goal will be to show that \Q m \ < \Pn\ on the contour, which by Rouche's Theorem gives us 
that P n and S have the same number of zeros inside the contour (which is none, since all of the roots of 
P n are positive real numbers). By taking contours bounding larger and larger regions contained in and 
tending towards the entire left-half plane of C, we will deduce that S has no roots with negative real 
component. 

We let our contour C — Cr consist of two parts: 

• a large half circle in the left-half of the complex plane 

• the portion of the imaginary axis connecting the two intersections of the above circle with the 
positive imaginary axis. 

First, we verify that \Q m \ < \P n \ on the half circle {|A| = R} in the left-half plane. Since P n is of 
higher order than Q m , as R — > oo, \P n \ dominates \Q m \. 

Next, we verify that \Q m \ < \P n \ on the imaginary axis. Again substituting A = iy for y > into 



(36), we have 



1^1 = ^J{y 2 +Pl)...{y 2 +Pl) (37) 

\Q m \ = biJ(y2 + q f)...( y 2 + q l). (38) 



Then, 



rfd (\Pn?\ I\U(V 2 +Pt) 



dy{\Q m \*) nr=i(y 2 +^ 2 )' [ ' 



and it follows that 



h ,± ( \p n ? \ = {nuy' + ^fcujiyiiu^+pt) 

dy\\Q m ?) (n™i(y 2 + ^ 2 )) 2 

(UtAy 2 +pf)) (S7=i (^ntwj^ + j) 



(40) 



Rewriting this, we find 

,o d np,\ 2 \ l 



dy\\Q m ?J UUv' + ih v=lV &=1 , Mi 




E \*y II (^+rf) -E br? n^+^ 



(41) 
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Ignoring the common denominator in both terms and considering only the first m terms of the first 
series, we obtain 

ml n \ rri / / 2 i 2\ n \ 

Eh* n m-e Q n &+&)• («> 

i=l y k=l,kfr ) j = l \\ y y i / i=l,i=£j J 

Now, comparing the terms of each series in order, we observe precisely the same pattern as in the 



main text; since pj < qj for all j = l...m, (42) is always greater than zero. Then, \P n \ > \Q m \ on 
the positive imaginary axis, and by symmetry also on the negative imaginary axis. Applying Rouchc's 
Theorem once again, we find that the number of roots of S — P n + Q m in the positive (negative) real 
half of the complex plane is identical to the number of roots of P n in the positive (negative) real half of 
the complex plane. 
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Figure Legends 




Figure 1. The semicirular contour Cr (dashed line) and the region it encloses (grey region) used in 
the proof of the cycle stability conjecture. By allowing the radius R of the semicircle to go to infinity, 
we are able to encompass the entire left-half of the complex plane. 
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Figure 2. Region of validity (in grey) for Stubborn Complex Roots Theorem. To apply the Stubborn 
Complex Roots Theorem, the complex roots of P n must have a smaller imaginary component than real 
component. 



